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ob ; ABSTRACT 

Nakamura, Sasaki, Tanaka, and Thorne have recently estimated the initial 
distribution of binary MACHOs in the galactic halo assuming that the MACHOs are 
I primordial black holes of mass ~ 0.5Mq, and considered their coalescence as a possible 

I source for ground- based interferometer gravitational wave detectors such as LIGO. 

Evolving their binary distribution forward in time to the present, the low- frequency 
J> ■ (lO^^Hz < / < lO^^Hz) spectrum of gravitational waves associated with such a 

population of compact binaries is calculated. The resulting gravitational waves would 
. form a strong stochastic background in proposed space interferometers such as LISA 

. and OMEGA. Low frequency gravitational waves are likely to become a key tool 

00 ! for determining the properties of binaries within the dark MACHO population. ( 

submitted to ApJ Letters September 3, 1998) 
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bJQ. Subject headings: gravitation — black holes — dark matter - — gravitational lensing - 



Galaxy: halo 



1. Introduction 

The MACHO project has thus far detected 15 candidate microlensing events in the direction 
of the Large Magellanic Cloud (LMC) ( |Cook et al. 1998 ) and 2 in the direction of the Small 



Magellanic Cloud (SMC) ( |Alcock et al. 1998| ). One of the LMC events ( [Bennett et al. 1996| ) and 



one of the SMC events ( [Alcock et al. 19"98| ) appear to have been caused by binary MACHO lenses. 
The LMC observations suggest that perhaps half of the halo consists of MACHOs of average mass 
~ 0.5Mq. Observational constraints make it unlikely that these MACHOs are main sequence 
red dwarfs ( [Bahcall et al. 1994 ; Graff fc Freese 1996a| ; Graff Sz Freese 1996b| ), or white dwarfs 



( Chariot fc Silk 1995 ; Adams fc Laughlin 1996| ). Given these difficulties, more speculative forms 



of MACHOs have been seriously considered, such as primordial black holes or boson stars. 

Recently, Nakamura, Sasaki, Tanaka, and Thorne (hereafter, NSTT) considered the possibility 
that the halo MACHOs consist of primordial black holes (BHMACHOs) ( [Nakamura et al. 1997| ). 



They derived an estimate of the initial binary distribution function, and determined the rate 
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at which BHMACHOs in binaries would coalesce due to emission of gravitational radiation. 
They found extragalactic halo BHMACHO coalescences to be promising hypothetical sources of 
gravitational waves (GW) for ground based interferometer detectors such as LIGO. 

In this letter BHMACHOs are again examined as possible sources of gravitational radiation, 
but now in the low frequency (10~^ — 10~^Hz) band where space-based interferometers, such as 
the proposed LISA ( Pender et al. 1998| ) and OMEGA ( Pelhngs et al. 1998| ) instruments wih be 
most sensitive. The NSTT binary distribution function is assumed, and evolved forward from the 
time of black hole formation to the present, to determine the distribution of BHMACHO binaries 
in semimajor axis and eccentricity today. The black hole binaries with orbital frequencies higher 
than about 10^^'^ Hz are in highly eccentric orbits, which greatly increases their efficiency as 
sources of gravitational radiation, and spreads that radiation over a large number of harmonics of 



the orbital frequency (Peters & Mathews 1963). Binary BHMACHOs are shown to create a strong 
stochastic background in the low frequency GW band, possibly dominating over all other galactic 
background sources. 



2. Evolution of the NSTT binary distribution 

The distribution of black hole MACHO binaries obtained by NSTT defines the probability 
distribution of values for the semimajor axis a and eccentricity e of the binary's orbit. Utilizing a 
simple model of primordial black hole formation, they found the probability distribution to be 

F{a,e)dade = {3/2)a^/'^x-^/^e{l - e^y^^^ dade, (1) 
where x = 1.1 x 10^*^ (Mb// /M©)^/^ (l^/i^) -4/3 

cm, and, as usual, 17 is the cosmological density 
parameter and h is the Hubble constant in units of lOOkms"^ Mpc~^. Allowable values for the 
semimajor axis are bounded above by x, and the eccentricity for each value of a is also bounded 
above, by Cmax = [1 - {a/xf/^]'^/^. 

The spatial distribution of black hole MACHO binaries is assumed to follow the standard 
spherical flat rotation halo model 

P = P~, 2 ' 

J.Z _|_ 

where p = O.OO79M0 pc~^ is the local density of dark matter, r is the Galactocentric radius, 
R = 8.5 kpc is the Galactocentric radius of the Sun, and a = 5.0 kpc is the halo core radius. 
The distribution of Eq.(|2|) is valid from r = out to some halo radius, Rh', typical values are 
Rh = 50 kpc, for a halo extending to the LMC; or R^ = 300 kpc for a halo extending half way to 
M31. 

The NSTT distribution given by Eq.(|l|) is very strongly peaked towards binaries with large 
eccentricities, due to the (1 — e^)^^/^ factor in F. For example, for binaries with initial orbital 
frequencies of ~ 10"^'^ Hz (which will evolve into the low frequency band today) 90% of the 
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binaries are within 5e ~ 0.00198 of the maximum eccentricity, Cmax — 0.9999775. In view of this, 
a simphfying assumption will be made: that all black hole MACHO binaries initially have the 
e = Cmax- The distribution of Eq.(|^) may then be rewritten in simplified form, after an integration 
over e, as: 

3 



F{a)da 2 



a\3/4 /a\3/2 



X J \x 



- (3) 
a 



The approximated NSTT distribution (with all binaries at their maximum value of e) 
describes the binary black hole MACHO population at the time of its formation, in the early 
universe. In order to calculate the spectrum of gravitational radiation produced by the binary 
black holes, it is necessary to know the distribution of binaries today. The population will 
have evolved due to a variety of effects, such as accretion, gravitational interactions with other 
MACHOs and with ordinary matter, and due to the emission of gravitational radiation. Only this 
last effect will be considered here. 

In order to determine the probability distribution of black hole MACHO binary orbital 
parameters today, the initial distribution given in Eq.(|3|) must be evolved forward in time from 
the early universe to the present epoch. The evolution of a binary system of point masses due 



to gravitational radiation reaction in the weak field limit is well understood ( Peters &: Mathews 



1963 ; Peters 1964 ). The distribution was evolved by numerically integration using a fourth-order 
adaptive stepsize Runge-Kutta routine. The results obtained are characterized by two functions: 
ao(aj), the present semimajor axis as a function of the initial semimajor axis (a subscript 
indicates the present day value of a quantity); and eo(ai), the present eccentricity as a function of 
the initial semimajor axis. 

The form of the distribution today is illustrated in Figure (1), in which the present eccentricity, 
Co, is plotted against the present orbital frequency, /q. The circularizing effect of gravitational 
radiation reaction is clearly visible; those binaries with orbital frequencies today of 10~^ Hz or 
below are still highly eccentric, while those which have evolved to higher frequencies such as 10"^ 
Hz are basically circular today. Integration of a 50 kpc halo with Qh? = 0.1 gives 2.6 x 10^ binaries 
in the low frequency band today, comparable to the expected number of close white dwarf binaries 
in the Galaxy (|Hils et al. 1990| ). In contrast, a 300 kpc halo with Q.h? = 1.0 will have 1.6 x 10^ 



black hole binaries in the low frequency band today. 



3. Stochastic gravitational wave background 

The calculation of the gravitational radiation spectrum emitted by the halo population of 
BHMACHO binaries differs substantially from the more familiar calculations for galactic white 
dwarf binaries. This is because the BHMACHO binaries will generally have non- negligible 
eccentricity. A binary system with an eccentric orbit will emit significantly more power in GW 
than a circular system with the same semimajor axis. Eccentric binaries will also emit GW on a 
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multitude of harmonics of the orbital frequency, whereas circular binaries emit purely in the n = 2 
mode. Peters and Mathews ( Peters fc Mathews 1963| ) found that the orbit-averaged power emitted 
in gravitational radiation by a binary system in the nth harmonic of the orbital frequency is 

32 ml ml (mi + ms) 

L{f,n,e) = - g{n,e) , 4 

5 

where / is the orbital frequency associated with semimajor axis a (not the frequency of the 
harmonic radiation), and g{n,e) is given by 



n 



q(n, e) = — 



Jn-2{ne) - 2eJn-i{ne) + -J„(ne) + 2eJ„+i(ne) - Jn+2{ne) 



n 



+ 1-eJ 



[Jn-2{ne) - 2Jn{ne) + Jn+2{ne)f + [J„(ne)]^| 



(5) 



The total GW signal at frequency / will incorporate contributions from a large number of 
binaries, and hence the signal will be stochastic in nature. The characteristic spectral amplitude 
at frequency / may be written as a sum over the uncorrelated contributions from the binaries for 
which / is a harmonic. 



1/2 
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where {dr"^) is the inverse distance squared to the source, averaged over the source distribution, 
and the frequency fj is the orbital frequency of a binary system for which the jth harmonic is at 
frequency: fj = f / j. Similarly, Cj is the eccentricity of the NSTT evolved binary for which the jth 
harmonic is at frequency /. The upper limit on the summation, jmax; is the value of j for which 
the semimajor axis of the binary is at the maximum value ~ x. For each value of j, dN/dfj 
gives the number of sources per unit frequency interval at frequency fj, and L{fj,j) is the the 
power emitted at frequency / by a binary with orbital frequency fj. 

The average squared inverse distance between halo sources and the Sun may be found 
by integrating over the source distribution given in Eq. (|^). For a 50 kpc halo, this gives 
((i-2)-V2 = 16.8 kpc, while for a 300 kpc halo (d^2^"i/2 ^ gg g 

The procedure followed to determine the low frequency GW spectrum of the NSTT 
distribution of black hole MACHO binaries is as follows. First, a frequency / is chosen. This fixes 
the frequencies (and hence semimajor axes) of the binaries whose jth harmonic will contribute to 
the sum in Eq.(^. For each harmonic index j in Eq.(|6|), the luminosity L(fj,j) and the number 
of binaries per unit frequency interval, dN/dfj, must then be determined. 

In order to determine the luminosity from Eq.(^), the harmonic index j must be known, 
along with the semimajor axis, a, and the eccentricity e. The semimajor axis is trivially obtained 
from the orbital frequency fj, while the eccentricity is also a function of fj, determined by the 
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numerical evolution of the approximate NSTT, as illustrated in Fig.(|T|). Once the eccentricity is 
known, the function g{j, e) may be evaluated. Direct numerical computation of Bessel function 
values was used for the lowest ten harmonics (j < 10). Accurate evaluation of the Bessel functions 
becomes more difficult for larger values of j. An asymptotic expansion was used to fix the value of 
g{n, 1) (which depends only on Jn{n)) and the shape of g{n,e) for arbitrary e was then modeled 
analytically using a combination of exponential and trigonometric functions. Trial variation of 
parameters and functional forms within the analytic fit indicate that the resulting luminosity is 
robustly estimated to about one percent accuracy. 

The evaluation of the number of binaries present today per unit frequency interval, dN/df, 
for an arbitrary binary frequency /, may be split into three pieces: 

dN _ ^ dF^ dai_ doo 
df dtti dao df 

Here dF/dai is the initial fraction of all black hole MACHO binaries per unit semimajor axis 
interval, dui/dao is the amount by which an infinitesimal interval in semimajor axis has expanded 
or contracted in evolving to the present day, and dao/df is the final conversion from semimajor axis 
to present orbital frequency via Kepler's third law. Given an orbital frequency /o today, the initial 
value of the semimajor axis Cj may be determined from the results of the numerical evolution 
of the approximate NSTT distribution. Then, dF/dui may be directly evaluated from Eq.(^). 
The evolution of dui/dao is also determined from the numerically evolved NSTT distribution. 
Multiplication of the factors in Eq.(|^ gives the number of binaries per unit frequency interval 
as a function of frequency, today. The results confirm the notion that the total GW signal from 
the black hole MACHO binaries will be stochastic in nature. Assuming an integration time of 
4 months (typical of LISA or OMEGA), every 10~^Hz bin will include modes from multiple 
binaries for all frequencies /^lO^^Hz. Above that frequency, individual bins may be occupied 
by individual binaries evolving towards coalescence, or may be empty, opening windows in which 
weaker sources (e.g., extragalactic) could be observed. 

The total gravitational wave spectral amplitude was then evaluated for a set of frequencies, /, 
incorporating contributions from harmonics up to a maximum value jmax- Setting jmax equal to 
1000 was found to provide an accuracy over the entire frequency range 10"^ - 10"^ Hz of better 
than one percent; specifically, the inclusion of 10^ harmonics rather than 10^ changes the spectral 
amplitude by less than 1% over this range. 

The resulting spectrum of gravitational radiation is shown in Figure (§). Both curves 
illustrated represent small halos, = 50 kpc, with O/i^ = 0.1 for the lower curve and = 1.0 
for the upper curve. The spectrum is only weakly dependent of the value of the halo radius; a 
larger halo simply adds additional sources at large distance. Curves for a 300 kpc halo are of 
the same shape as those illustrated, but with spectral amplitude increased by a constant factor 
of ~ 1.106 over the full frequency range (on the logarithmic plot of Fig(Q), this corresponds to 
moving the curves upwards by 0.0436). 
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The peak in the spectrum at around / = 10~^ Hz is due to a combination of circumstances. 
Binaries whose n = 2 mode is at that frequency have eccentricities of around 0.4 (see Fig.([l|)), 
so that still about 50% of their power is being emitted in the n = 2 mode ( Peters fc Mathews] 



1963). At the same time, lower frequency binaries within an order of magnitude in frequency are 
significantly eccentric, so that they contribute substantial power at 10^'^ Hz through their higher 
harmonics. At higher frequencies, the lower frequency binaries are more nearly circular, and 
their higher harmonics do not contribute so much power. At lower frequencies, the fundamental 
mode's contribution (and, later, the higher modes' contributions) drops off sharply due to the 
dependence in Eq.(Q). 

At high frequencies, the spectrum is seen to approach a power law behavior, /i/ ~ f~'^^^- 
Since binaries in the higher frequency range are nearly circular, essentially all of the signal 
is generated by the n = 2 harmonic. The observed power law is then explained by noting 
that L ~ f^'^^^, while the factors of dN/df in Eq.(^ have the following asymptotic form: 
dF/dai ~ constant, da j/doo ~ f~'^,dao/df ~ f^^^^- Using these asymptotic forms in Eq.( ^) and 
only including the j = 2 mode yields /i/ ~ /~^^^- 

Also shown in Fig.(0) are the sensitivity curves for the proposed LISA (Bender et al. 1998) 



and OMEGA ( Pellings et al. 1998 ) space interferometer GW detectors, and a recent estimate 



of the background due to galactic close white dwarf binaries (CWDB) (|Hils fc Bender 1998 ) 



Except at low frequencies, the BHMACHO binary background would be notably stronger than the 
CWDB background. The halo black hole binaries would form a strong "confusion noise" stochastic 
background, significantly larger than the instrument noise in either proposed interferometer over a 
large frequency range 10~^Hz < / < lO^^Hz. Of course, while the stochastic signal is in one sense 
a noise, it also provides direct information on the nature and distribution of MACHO binaries in 
the Galactic halo. 

Whether the halo MACHOs are primordial black holes, white dwarfs, boson stars, or some 
other type of compact object, some fraction of them will be in binary systems. Given the huge 
number of MACHOs in the Galaxy's halo, the low frequency gravitational wave signal from 
these binary systems will quite likely be detectable by the space-based interferometers. These 
interferometers will be excellent instruments for determining the nature and population of the 
halo MACHO binaries. 
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Fig. 1. — The eccentricity eo of the evolved NSTT black hole MACHO binaries is plotted against 
their orbital frequency /o today. 



Fig. 2. — The spectral amplitude of gravitational waves emitted by an evolved NSTT distribution of 
BHMACHO binaries in a 50 kpc Galactic halo is shown for ^Ih^ = 0.1 (lower curve) and Qh'^ = 1.0 
(upper curve). Also shown are the expected sensitivities of the proposed LISA and OMEGA space- 
based interferometers, and a recent estimate of the stochastic background due to Galactic close 
white dwarf binaries (curve labeled CWDB). 



